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ABSTRACT 



in ■ 
o 

f"*^ ' Context. Observations of very high energy y-rays from blazars provide information about acceleration mechanisms occurring in their innermost 
regions. Studies of variability in these objects allow a better understanding of the mechanisms at play. 
• • . Aims. To investigate the spectral and temporal variability of VHE (> 100 GeV) y-rays of the well-known high-frequency-peaked BLLac object 
^ PKS 2155-304 with the H.E.S.S. imaging atmospheric Cherenkov telescopes over a wide range of flux states. 

' Methods. Data collected from 2005 to 2007 are analyzed. Spectra are derived on time scales ranging from 3 years to 4 minutes. Light curve 
rS ] variability is studied through doubling timescales and structure functions, and is compared with red noise process simulations. 

Results. The source is found to be in a low state from 2005 to 2007, except for a set of exceptional flares which occurred in July 2006. The 
quiescent state of the source is characterized by an associated mean flux level of (4.32 ± 0.09 sta t ± 0.86 syst ) x 10~" cirT 2 s _1 above 200 GeV, or 
approximately 15% of the Crab Nebula, and a power law photon index of T = 3.53 ± 0.06 stat ± 0. 10 syst . During the flares of July 2006, doubling 
timescales of ~ 2 min are found. The spectral index variation is examined over two orders of magnitude in flux, yielding different behaviour at low 
and high fluxes, which is a new phenomenon in VHE y-ray emitting blazars. The variability amplitude characterized by the fractional r.m.s. F var is 
strongly energy-dependent and is oc £°- 19±0 - 01 . l^jgj^j curve r.m.s. correlates with the flux. This is the signature of a multiplicative process which 
can be accounted for as a red noise with a Fourier index of ~ 2. 

Conclusions. This unique data set shows evidence for a low level y-ray emission state from PKS 2155-304, which possibly has a different origin 
than the outbursts. The discovery of the light curve lognormal behaviour might be an indicator of the origin of aperiodic variability in blazars. 

Key words, gamma rays: observations - Galaxies : active - Galaxies : jets - BL Lacertae objects: individual objects: PKS 2155-304 
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getic and extremely variable extragalactic very high energy y-ray 
emitters. Their luminosity varies in unpredictable, highly irreg- 
ular ways, by orders of magnitude, and at all wavelengths across 
the electromagnetic spectrum. The very high energy (VHE, E > 
100 GeV) y-ray fluxes vary often on the shortest timescales that 
can be seen in this type of object, with large amplitudes which 
can dominate the overall output. It hence indicates that the un- 
derstanding of this energy domain is the most important one for 
understanding the underlying fundamental variability and emis- 
sion mechanisms at play in high flux states. 

It has been, however, difficult to ascertain whether y-ray 
emission is present only during high flux states or also when the 
source is in a more stable or quiescent state but with a flux which 
is below the instrumental limits. The advent of the current gener- 
ation of atmospheric Cherenkov telescopes with unprecedented 
sensitivity in the VHE regime gives new insights into these ques- 
tions. 

The high frequency peaked BL Lac object (HBL) 
PKS 2155-304, located at a redshift z = 0.117, initially dis- 
covered as a VHE y-ray emitter by the Mark 6 telescope 
dChadwick et al. 1999k h as been detected by the first H.E.S.S. 
telescope in 2002-2003 (Aha ronian et al.l 2005b). It has been 
frequently observed by the full array of four telescopes since 
2004, either sparsely during the H.E.S.S. monitoring program, 
or intensely during dedicated campaigns such as that de- 
scribed in lAharonian et al.l (2005c), showing mean flux lev- 
els of ~ 20% of the Crab Nebula flux for energies above 
200 GeV. During the summer of 2006, PKS 2155-304 exhib- 
ited unprecedented flux levels accompanied by strong variability 
(lAharonian et al. 1 2007a), making temporal and spectral variabil- 
ity studies possible on timescales of the order of a few minutes. 
The VHE y-ray emission is usually thought to originate from a 
relativistic jet, emanating from the vicinity of a Supermassive 
Black Hole (SMBH). The physical processes at play are still 
poorly understood, but the analysis of the y-ray flux spectral and 
temporal characteristics is well suited to provide better insights. 

For this goal the data set of H.E.S.S. observations of 
PKS 2155-304 between 2005 and 2007 is used. After describing 
the observations and the analysis chain in Section|2] the emission 
from the "quiescent", i.e. nonflaring, state of the source will be 
characterized in Section [3] Section |4] details spectral variability 
related to the source intensity. Section will focus on the de- 
scription of the temporal variability during the highly active state 
of the source, and its possible energy dependence. Section[6]will 
illustrate a description of the observed variability phenomenon 
by a random stationary process, characterized by a simple power 
density spectrum. Section|7]will show how limits on the charac- 
teristic time of the source can be derived. The multi-wavelength 
aspects from the high flux state will be presented in a second 
paper. 

2. Observations and analysis 

H. E.S.S. is an array of four imaging atmospheric Cherenkov 
telescopes situated in the Khomas Highland of Namibia 
(23°16'18" South, 16°30'00" East ), at an elevat ion of 

I, 800 meters above sea level (see Aharoni an et al.l 2006). 
PKS 2155-304 was observed by H.E.S.S. each year since 
2002; results of observations in 2002, 2003 and 2004 can be 
found in lAharonian et al.l (2005b). lAharonian et al.l (2005c) and 
|Giebels et al. (200"5j| The data reported here were collected be- 
tween 2005 and 2007. In 2005, 12.2 hours of observations were 
taken. A similar observation time was scheduled in 2006, but 
following the strong flare of July 26th ( lAharonian et al.l 2007a) 



Table 1. Summary of observations for each year. T represents 
the live-time (hours), « on the number of on-source events, n s 
the number of off-source events (from a region five times larger 
than for the on-source events), and <x the significance of the cor- 
responding excess, given in units of standard deviations. 



Year 


T 




"off 


Excess 


a 


cr/VT 


2005 


9.4 


7,282 


27,071 


1,868 


21.8 


7.1 


2006 


66.1 


123,567 


203,815 


82,804 


288.4 


35.5 


2007 


13.8 


11,012 


40,065 


2,999 


28.6 


7.7 


Total 


89.2 


141,861 


270,951 


87,671 


275.6 


29.2 




zenith angle (cleg.) 

Fig. 1. Zenith angle distribution for the 202 4-telescopes ob- 
servation runs from 2005 to 2007. The inset shows, for each 
zenith angle, the energy threshold associated with the analysis 
presented in Section[2] 



it was decided to increase this observation time significantly. 
Ultimately, from June to October 2006, this source was observed 
for 75.9 hours, with a further 20.9 hours in 2007. 

The data have been recorded during runs of 28 minutes nom- 
inal duration, with the telescopes pointing at 0.5° from the source 
position in the sky to enable a simultaneous estimate of the back- 
ground. This offset has been taken alternatively in both right as- 
cension and declination (with both signs), in order to minimize 
systematics. Only the runs passing the H.E.S.S. data-quality se- 
lection criteria have been used for the analyses presented below. 
These criteria imply good atmospheric conditions and checks 
that the hardware state of the cameras is satisfactory. The num- 
ber of runs thus selected is 22 for 2005, 153 for 2006 and 35 
for 2007, corresponding to live-times of 9.4, 66.1 and 13.8 hours 
respectively. During these observations, zenith angles were be- 
tween 7 and 60 degrees, resulting in large variations in the in- 
strument energy threshold (E±, see Fig. [1} and sensitivity. This 
variation has been accounted for in the spectral and temporal 
variability studies presented below. 

The data have been analysed following the prescription pre- 
sented in Ah aronian et al.l (2006), using the loose set of cuts 
which are well adapted for bright sources with moderately soft 
spectra, and the Reflected-Region method for the definition of the 
on-source and off-source data regions. A year-wise summary of 
the observations and the resulting detections are shown in Table 
Q] A similar summary is given in AppendixlAlfor the 67 nights of 
data taken, showing that the emission of PKS 2155-304 is easily 
detected by H.E.S.S. almost every night. For 66 nights out of 67, 
the significance per square root of the live-time (tx/ yff, where 
T is the observation live-time) is at least equal to 3.6 cr/ ^fh, the 
only night with a lower value — MJD 53705 — corresponding 
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to a very short exposure. In addition, for 61 nights out of 67 
the source emission is high enough to enable a detection of the 
source with 5 cr significance in one hour or less, a level usually 
required in this domain to firmly claim a new source detection. In 
2006 the source exhibits very strong activity (38 nights, between 
MJD 53916 - 53999) with a nightly cr/ Vr varying from 3.6 to 
150, and being higher than 10 cr/ yfh for 19 nights. The activity 
of the source climaxes on MJD 53944 and 53946 with statistical 
significances which are unprecedented at these energies, the rate 
of detected y-rays corresponding to 2.5 and 1.3 Hz, with 150 
and 98 cr/ Vh respectively. 

For subsequent spectral analysis, an improved energy re- 
construction method with respect to the one described in 
lAharonian et al.l (2006) was applied. This method is based on a 
look-up table determined from Monte-Carlo simulations, which 
contains the relation between an image's amplitude and its re- 
constructed impact parameter as a function of the true energy, 
the observation zenith angle, the position of the source in the 
camera, the optical efficiency of the telescopes (which tend to 
decrease due to the aging of the optical surfaces), the number of 
triggered telescopes and the reconstructed altitude of the shower 
maximum. Thus, for a given event, the reconstructed energy is 
determined by requiring the minimal x 1 between the image am- 
plitudes and those expected from the look-up table correspond- 
ing to the same observation conditions. This method yields a 
slightly lower energy threshold (shown in Fig. Q] as a function 
of zenith angle), an energy resolution which varies from 15% to 
20% over all the energy range, and biases in the energy recon- 
struction which are smaller than 5%, even close to the threshold. 
The systematic uncertainty in the normalization of the H.E.S.S. 
energy scale is estimated to be as large as 15%, corresponding 
for such soft spectrum source to 40% in the overall flux normal- 
ization as quoted in Aha ronian et al.l (2009). 

All the spectra presented in this paper have been obtained 
using a forward-folding maximum likelihood method based on 
the measured energy-dependent on-source and off-source distri- 
butions. This method, fully described in |Piron et al. (200"T)] per- 
forms a global deconvolution of the instrument functions (energy 
resolution, collection area) and the parametrization of the spec- 
tral shape. Two different sets of parameters, corresponding to a 
power law and to a power law with an exponential cut-off, are 
used for the spectral shape, with the following equations : 

<ke) = <po{§^y r (i) 

0(£) = 0oexp(-^)(-^)~ r exp(--^) (2) 

^-cut t^Q ^cut 

4>o represents the differential flux at Eq (chosen to be 1 TeV), T is 
the power law index and E cut the characteristic energy of the ex- 
ponential cut-off. The maximum likelihood method provides the 
best set of parameters corresponding to the selected hypothesis, 
and the corresponding error matrix. 

Finally, various data sets have been used for subsequent anal- 
yses. These are summarized in Table|2] 

3. Characterization of the quiescent state 

As can be seen in Fig. [2] with the exception of the high state 
of July 2006 PKS 2155-304 was in a low state during the ob- 
servations from 2005 to 2007. This section explores the vari- 
ability of the source during these periods of low-level activity, 
based on the determination of the run-wise integral fluxes for the 
data set Dqs , which excludes the flaring period of July 2006 and 
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Fig. 2. Monthly averaged integral flux of PKS 2155-304 above 
200 GeV obtained from data set D (see Table [2}, which corre- 
sponds to the 165 4-telescope runs whose energy threshold is 
below 200 GeV. The dotted line corresponds to 15% of the Crab 
Nebula emission level (see Section [3~!2b . 

also those runs whose energy threshold is higher than 200 GeV 
(see l3.1l for justification). As for Sections[5]and [6] the control of 
systematics in such a study is particularly important, especially 
because of the strong variations of the energy threshold through- 
out the observations. 

3. 1 . Method and systematics 

The integral flux for a given period of observations is determined 
in a standard way. For subsequent discussion purposes, the for- 
mula applied is given here : 

J^S(E)dE 

* = r oo rE ( 3 ) 

T J Q J E m " A(E)R(E,E')S(E)dE'dE 

where T represents the corresponding live-time, A(E) and 
R(E, E') are, respectively, the collection area at the true energy 
E and the energy resolution function between E and the mea- 
sured energy E' , and S(E) the shape of the differential energy 
spectrum as defined in Eq.Q]and[2]. Finally, N exp is the number 
of measured events in the energy range [E m j n ,E max ]. 

In the case that S (E) is a power law, an important source of 
systematic error in the determination of the integral flux vari- 
ation with time comes from the value chosen for the index F. 
The average 2005-2007 energy spectrum yields a very well de- 
termined power law indejfl However, in Section [4] it will be 
shown that this index varies depending on the flux level of the 
source. Moreover, in some cases the energy spectrum of the 
source shows some curvature in the TeV region, giving slight 
variations in the fitted power law index depending on the energy 
range used. 

For runs whose energy threshold is lower than E m [ n , a sim- 
ulation performed under the observation conditions correspond- 
ing to the data shows that an index variation of Ar = 0.1 implies 
a flux error at the level of AO ~ 1 %, this relation being quite lin- 
ear up to AF ~ 0.5. However, this relation no longer holds when 

1 The average 2005-2007 energy spectrum yields a power law with a 
photon index T = 3.37 ± 0.02 stat . One should note that some curvature is 
observed at higher energies, resulting in a better spectral determination 
when the alternative hypothesis shown in (Eq. [2]l is used, yielding a 
harder index (r = 3.05 ± 0.05 stat ) with an exponential cut-off at energy 
Ecut = 1.76±0.27 sta t TeV. This curved model is prefered at a 8.4 cr level 
as compared to the power law hypothesis. However, the choice of the 
model has little effect on the determination of the integral flux values 
above 200 GeV, the integral being dominated by the low-energy part of 
the energy spectrum. 
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Table 2. The various data sets used in the paper, referred to in the text by the labels presented in this table. Only runs with the 
full array of four telescopes in operation (202 runs over 210) and an energy threshold lower than 200 GeV (165 runs over 202) 
are considered. The corresponding period of the observations, the number of runs, the live-time T (hours), the number of y excess 
events and its significance cr are shown. The column section indicates the sections of the paper in which each data set is discussed. 
Additional criteria for the data set definitions are indicated in the last column. 



Label 


Period 


Runs 


T(hours) 


Excess 


cr 


Section 


Additional criteria 


D 


2005-2007 


165 


69.7 


67,654 


237.4 


|4J|7| 




D QS 


2005-2007 


115 


48.1 


12,287 


60.5 


133113^1 1331171 


July 2006 excluded 


Dqs -2005 


2005 


19 


8.0 


1,816 


22.6 


m 




^gS-2006 


2006 


61 


26.3 


7,472 


48.4 




July 2006 excluded 


^CW-2007 


2007 


35 


13.8 


2,999 


28.6 


E3 




Djuly06 


July 2006 


50 


21.6 


55,367 


281.8 






Dflares 


July 2006 (4 nights) 


27 


11.8 


46,036 


284.1 


SI3|6l|7j 





the energy threshold is above E m [ n , as the determination of <1> 
becomes much more dependent on the choice of F. For this rea- 
son, only runs whose energy threshold is lower than E m [ n will be 
kept for the subsequent light curves. The value of E m [ n is chosen 
as 200 GeV, which is a compromise between a low value which 
maximises the excess numbers used for the flux determinations 
and a high value which maximises the number of runs whose 
energy threshold is lower than E m \ n . 

3.2. Run-wise distribution of the integral flux 

From 2005 to 2007, PKS 2155-304 is almost always detected 
when observed (except for two nights for which the exposure 
was very low), indicating the existence, at least during these ob- 
servations, of a minimal level of activity of the source. Focussing 
on data set Dqs (which excludes the July 2006 data where the 
source is in a high state), the distribution of the integral fluxes 
of the individual runs above 200 GeV has been determined for 
the 1 15 runs, using a spectral index T = 3.53 (the best value for 
this data set, as shown in 13.4b . This distribution has an asym- 
metric shape, with mean value (4.32 + 0.09 sta t) X 10~ n cm _2 s _1 
and root mean square (r.m.s.) (2.48 + 0.11 st at) x 10~ n cm _2 s _1 , 
and is very well described with a lognormal function. Such a 
behaviour implies that the logarithm of fluxes follows a nor- 
mal distribution, centered on the logarithm of (3.75 + 0.11 sta t) x 
10" cm 2 s~'. This is shown in the left panel of Fig. [3] where 
the solid line represents the best fit obtained with a maximum- 
likelihood method, yielding results independent of the choice 
of the intervals in the histogram. It is interesting to note that 
this result can be compared to the fluxes measured by H.E.S.S. 
from PKS 2155-304 during its construction phase, in 2002 and 
2003 (see lAharonian et al.l 2005b and lAharonian et al.l 2005c). 
As shown in Table [3] these flux levels extrapolated down to 
200 GeV were close to the value corresponding to the peak 
shown in the left panel of Fig. [3] 

The right panel of Fig. [3] shows how the flux distribution is 
modified when the July 2006 data are taken into account (data set 
D in Table [2}: the histogram can be accounted for by the super- 
position of two Gaussian distributions (solid curve). The results, 
summarized in Table [H are also independent of the choice of the 
intervals in the histogram. Remarkably enough, the character- 
istics of the first Gaussian obtained in the first step (left panel) 
remain quite stable in the double Gaussian fit. 

This leads to two conclusions. First, the flux distribution of 
PKS 2155-304 is well described considering a low state and a 
high state, for each of which the distribution of the logarithms 
of the fluxes follows a Gaussian distribution. The characteristics 
of the lognormal flux distribution for the high state are given 



Table 3. Integral fluxes and their statistical errors from 
2002 and 2003 observations of PKS 2155-304 during the 
H.E.S.S. construction phase. These values are taken from 
lAharonian et al.l (2005b) and lAharonian et al.l (2005c) and cor- 
respond to flux extrapolations to above 200 GeV. 



Month 


Year 


O [10-" 


cm 2 s '] 


July 


2002 


16.4 


±4.7 


Oct. 


2002 


8.9 


±5.2 


June 


2003 


5.8 


± 1.4 


July 


2003 


2.9 


±0.5 


Aug. 


2003 


3.5 


±0.5 


Sep. 


2003 


4.9 


± 1.2 


Oct. 


2003 


5.2 


±0.5 



Table 4. The distribution of the flux logarithm. First column: 
distribution as fitted by a single Gaussian law for the "quiescent" 
regime (data set Dqs ). Second column: distribution fitted by 
two Gaussian laws, one for the "quiescent" regime, the other for 
the flaring regime (data set D). Decimal logarithms are quoted 
to make the comparison with the left panel of Fig. [3]easier and 
the flux is expressed in cirr 2 s . In the first line the average of 
fluxes is reported, while in the second line their r.m.s.. 



"Quiescent" regime Flaring regime 
(log 10 O) -10.42 ± 0.02 -9.79 ±0.11 
r.m.s. of logip O 0.24 ± 0.02 0.58 ± 0.04 



in Sections |5] and |7J Secondly, PKS 2155-304 has a level 
of minimal activity which seems to be stable on a several-year 
time-scale. This state will henceforth be referred to as the 
"quiescent state" of the source. 



3.3. Width of the run-wise flux distribution 

In order to determine if the measured width of the flux distribu- 
tion (left panel of Fig. can be explained as statistical fluctua- 
tions from the measurement process a simulation has been car- 
ried out considering a source which emits an integral flux above 
200 GeV of 4.32 x 10~ n cm _2 s _1 with a power law photon spec- 
trum index F = 3.53 (as determined in the next section). For 
each run of the data set Dqs the number n y expected by convolv- 
ing the assumed differential energy spectrum with the instrument 
response corresponding to the observation conditions is deter- 
mined. A random smearing around this value allows statistical 



HESS Collaboration et al.: VHE y-ray emission of PKS 2155-304: spectral and temporal variability 



5 




log 10 <E> (cm s 1 ) log 10 * (cm s 1 ) 

Fig. 3. Distributions of the logarithms of integral fluxes above 200 GeV in individual runs. Left: from 2005 to 2007 except the July 
2006 period (data set Dqs), fitted by a Gaussian. Right: all runs from 2005 to 2007 (data set D), where the solid line represents the 
result of a fit by the sum of 2 Gaussians (dashed lines). See Table[4]fbr details. 



fluctuations to be taken into account. The number of events in 
the off-source region and also the number of background events 
in the source region are derived from the measured values « ff in 
the data set. These are also smeared in order to take into account 
the expected statistical fluctuations. 

10000 such flux distributions have been simulated, and for 
each one its mean value and r.m.s. (which will be called be- 
low RMSD) are determined. The distribution of RMSD thus ob- 
tained, shown in Fig. [4] is well described by a Gaussian centred 
on 0.98 x 10~ n cm _2 s _1 (which represents a relative flux disper- 
sion of 23%) and with a (Trmsd of 0.07 x 10~ n cirrV 1 . 

It should be noted that here the effect of atmospheric fluc- 
tuations in the determination of the flux is only taken into ac- 
count at the level of the off-source events, as these numbers are 
taken from the measured data. But the effect of the correspond- 
ing level of fluctuations on the source signal is very difficult 
to determine. If a conservative value of 20% is considered Q, 
which is added in the simulations as a supplementary fluctua- 
tion factor for the number of events expected from the source, 
a RMSD distribution centred on 1.30 x 10~ n crrT^ 1 with a 
crmsd of 0.09 x 10~ n cm^s -1 is obtained. Even in this con- 
servative case, the measured value for the flux distribution r.m.s. 
((2.48 ± 0.1 l s tat) x 10~ n cirrV) is very far (more than 8 stan- 
dard deviations) from the simulated value. All these elements 
strongly suggest the existence of an intrinsic variability associ- 
ated with the quiescent state of PKS 2155-304. 

3.4. Quiescent-state energy spectrum 

The energy spectrum associated with the data set Dqs , shown in 
Fig. |5] is well described by a power law with a differential flux at 
lTeVof^o = (l.SliO.BstaJxlO^cirrV'TeV- 1 and an index 

2 A similar procedure has been carried out on the Crab Nebula ob- 
servations. Considering this source to be perfectly stable it allows us to 
determine an upper limit to the fluctuations of the Crab signal due to the 
atmosphere, and a value of ~ 10% was derived. Nonetheless, this value 
is linked to the observations' epoch and zenith angles, and to the source 
spectral shape. 




Fig. 4. Distribution of RMSD obtained when the instrument 
response to a fixed emission (<E> = 4.32 x 10~ n crrT 2 s~' and 
F = 3.53) is simulated 10000 times with the same observation 
conditions as for the 115 runs of the left part of Fig. [3] 



of r = 3.53 ± 0.06 sta t- The stability of these values for spectra 
measured separately for 2005, 2006 (excluding July) and 2007 
is presented in Table [5] The corresponding average integral flux 
is (4.23 ± 0.09 sta t) X 10™ 11 cm _2 s _1 , which is as expected in very 
good agreement with the mean value of the distribution shown 
in the left panel of Fig. [3] 

Bins above 2 TeV correspond to y-ray excesses lower than 
20 y and significances lower than 2 cr. Above 5 TeV excesses are 
even less significant (~ lcr or less) and 99% upper-limits are 
used. There is no improvement of the fit when a curvature is 
taken into account. 
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Fig. 5. Energy spectrum of the quiescent state for the period 
2005-2007. The green band correponds to the 68% confidence- 
level provided by the maximum likelihood method. Points are 
derived from the residuals in each energy bin, only for illustra- 
tion purposes. See Section [3~4l for further details. 

Table 5. Parametrization of the differential energy spectrum 
of the quiescent state of PKS 2155-304, determined in the 
energy range 0.2-10 TeV, first for the 2005-2007 period and 
also separately for the 2005, 2006 (excluding July) and 2007 
periods. Corresponding data sets are those of Table [2] 0o 
(10~ 12 cirrV'TeV- 1 ) is the differential flux at 1 TeV, T the pho- 



ton index and O(10 
Errors are statistical. 



cm s ) the integral flux above 0.2 TeV. 



Year 


Data set 




h 




r 




<D 


2005-2007 


D QS 


1.81 


±0.13 


3.53 


±0.06 


4.23 


±0.09 


2005 


Dqs -2005 


1.59 


+ 0.32 


3.56 


±0.16 


3.83 


±0.21 


2006 


OgS-2006 


1.87 


+ 0.18 


3.59 


±0.08 


4.65 


±0.13 


2007 


Dqs -2007 


1.84 


±0.24 


3.43 


±0.11 


3.78 


±0.16 



4. Spectral variability 

4.1. Variation of the spectral index for the whole data set 
2005-2007 

The spectral state of PKS 2155-304 has been monitored since 
2002. The first set of observations (lAharonian et al.l 2005b), 
from July 2002 to September 2003, shows an average energy 
spectrum well described by a power law with an index of F = 
3.32+0.06 st at, for an integral flux (extrapolated down to 200 GeV 
) of (4.39 + 0.40 sta t)x 10 -11 cm^s"'. No clear indication of spec- 
tral variability was seen. Consecutive observations in October 
and November 2003 (lAharonian et al.l 2005c) gave a similar 
value for the index, F = 3.37 + 0.07 stat , for a slightly higher flux 
of (5.22 + 0.54 stat ) x IO" 11 cnrV 1 Later, during H .E.S.S. ob- 
servations of the first (MJD 53944, lAharonian et all 2007a) and 
second (MJD 53946, lAharonian et al.l 2009) exceptional flares 
of July 2006, the source reached much higher average fluxes, 
corresponding to (1.72 + 0.05 stat ) X 10 _9 cm _2 s _1 and (1.24 



0.02 stat )xl0 9 cm 2 s '0 respectively. In the first case, no strong 
indications for spectral variability were found and the average 
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Fig. 6. Integral flux above 200 GeV measured each night dur- 
ing late July 2006 observations. The horizontal dashed line 
corresponds to the quiescent state emission level defined in 
Section [HI 



index F = 3.19 ± 0.02 stat was close to those associated with the 
2002 and 2003 observations. In the second case, clear evidence 
of spectral hardening with increasing flux was found. 

The observations of PKS 2155-304 presented in this paper 
also include the subsequent flares of 2006 and the data of 2005 
and 2007. Hence, the evolution of the spectral index is studied 
for the first time for a flux level varying over two orders of mag- 
nitude. This spectral study has been carried out over the fixed 
energy range 0.2-1 TeV in order to minimize both systematic 
effects due to the energy threshold variation and the effect of 
the curvature observed at high energy in the flaring states. The 
maximal energy has been chosen to be at the limit where the 
spectral curvature seen in high flux states begins to render the 
power law or exponential curvature hypotheses distinguishable. 
As flux levels observed in July 2006 are significantly higher than 
in the rest of the data set (see Fig. the flux-index behaviour 
is determined separately first for the July 2006 data set itself 
{DjuLYoe) and secondly for the 2005-2007 data excluding this 
data set (Dqs). 

On both data sets, the following method was applied. The 
integral flux was determined for each run assuming a power law 
shape with an index of F = 3.37 (average spectral index for the 
whole data set), and runs were sorted by increasing flux. The set 
of ordered runs was then divided into subsets containing at least 
an excess of 1500 y above 200 GeV and the energy spectrum of 
each subset was determinecfl 

The left panel of Fig. [7] shows the photon index versus in- 
tegral flux for data sets Dqs (grey crosses) and D julyq6 (black 
points). Corresponding numbers are summarized in AppendixlBl 
While a clear hardening is observed for integral fluxes above 
a few 10~ 10 cm _2 s _1 , a break in this behaviour is observed for 
lower fluxes. Indeed, for the data set Z)/j/lfo6 (black points) 
a linear fit yields a slope dr/dO = (3.0 ± 0.3 stat ) x 10 8 cm 2 s, 
whereas the same fit for data set Dqs (grey crosses) yields a 
slope dF/dd) = (-3.4 ± 1.9 stat ) x 10 9 cm 2 s. The latter corre- 
sponds to a^- 2 probability POf 2 ) = 71%; a fit to a constant yields 
P(x 2 ) = 33% but with a constant fitted index incompatible with 
a linear extrapolation from higher flux states at a 3 cr level. This 
is compatible with conclusions obtained either with an indepen- 
dant analysis or when these spectra are processed following a 
different prescription. In this prescription the runs were sorted 
as a function of time in contiguous subsets with similar photon 
statistics, rather than as a function of increasing flux. 



corresponding to data set T200 in Aharonian et al. 1 (2009) 



4 even for lower fluxes, the significance associated with each subset 
is always higher than 20 standard deviations 
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Fig. 7. Evolution of the photon index T with increasing flux O in the 0.2-1 TeV energy range. The left panel shows the results for 
the July 2006 data (black points, data set Djuiyoe) and for the 2005-2007 period excluding July 2006 (grey points, data set Dqs). 
The right panel shows the results for the four nights flaring period of July 2006 (black points, data set Deiares) an d one point 
corresponding to the quiescent state average spectrum (grey point, again data set Dqs). See text in Sections [4. II (left panel) and 14.21 
(right panel) for further details on the method. 



The form of the relation between the index versus inte- 
gral flux is unprecedented in the TeV regime. Prior to the 
results presented here, spectral variability has been detected 
only in two other blazars, Mrk421 and Mrk501. For Mrk421, 
a clear hardening with increasing flux appeared during the 
1999/2000 and 2000/20 01 observations performed with HEGRA 
dAharonian et al. 20021) and also during the 2004 observations 
performed with H.E.S.S. dAharonian e t al. 2005a). In addition, 
the Mrk501 observations carried out with CAT during the strong 
flares of 1997 ( |Djannati-Atai et al. 1999 1 and also the recent 
observation performed by MAGIC in 2005 ( I Albert et al.1 2007) 
have shown a similar hardening. In both studies, the VHE peak 
has been observed in the vF v distributions of the flaring states of 
Mrk501. 



4.2. Variation of the spectral index for the four flaring nights 
of July 2006 

In this section, the spectral variability during the flares of July 
2006 is described in more detail. A zoom on the variation of 
the integral flux (4-minute binning) for the four nights contain- 
ing the flares (nights MJD 53944, 53945, 53946, and 53947, 
called the "flaring period") is presented in the top panel of 
Fig. [8] This figure shows two exceptional peaks on MJD 53944 
and MJD 53946 which climax respectively at fluxes higher than 
2.5 x lO^cirrV 1 and 3.5 x lO^cnrV 1 (~ 9 and ~ 12 times 
the Crab Nebula level above the same energy), both about two 
orders of magnitude above the quiescent state level. 

The variation with time of the photon index is shown in 
the bottom panel of Fig. [8] To obtain these values, the y ex- 
cess above 200 GeV has been determined for each 4-minute bin. 
Then, successive bins have been grouped in order to reach a 
global excess higher than 600 y. Finally, the energy spectrum 
of each data set has been determined in the 0.2-1 TeV en- 
ergy range, as before (corresponding numbers are summarized 
in Appendix Table IB. 4b . There is no clear indication of spectral 
variability within each night, except for MJD 53946 as shown in 
lAharonian et al. 1 (2009). However, a variability can be seen from 
night to night, and the spectral hardening with increasing flux 



level already shown in Fig. [7] is also seen very clearly in this 
manner. 

It is certainly interesting to directly compare the spectral 
behaviour seen during the flaring period with the hardness of 
the energy spectrum associated with the quiescent state. This is 
shown in the right panel of Fig. [7j where black points corre- 
spond to the four flaring nights; these were determined in the 
same manner as for the left panel (see 14. II for details). A linear 
fit here yields a slope dr/d<D = (2.8 ±0.3 sta t) x 10 8 cm 2 s. The grey 
cross corresponds to the integral flux and the photon index asso- 
ciated with the quiescent state (derived in a consistent way in the 
energy range from 0.2-1 TeV), showing a clear rupture with the 
tendancy at higher fluxes (typically above 10~ 10 cm _2 s _1 ). 

These four nights were further examined to search for dif- 
ferences in the spectral behaviour between periods in which 
the source flux was clearly increasing and periods in which it 
was decreasing. For this, the first 16 minutes of the first flare 
(MJD 53944) are of special interest because they present a very 
symmetric situation: the flux increases during the first half, and 
then decreases to its initial level, the averaged fluxes are sim- 
ilar in both parts (~ 1.8 x 10 -9 cm _2 s _1 ), and the observation 
conditions (and thus the instrument response) are almost con- 
stant — the mean zenith angle of each part being respectively 
7.2 and 7.8 degrees. Again, the spectra have been determined in 
the 0.2-1 TeV energy range, giving indices of T — 3.27 + 0.12 stat 
and T = 3.27 + 0.09 sta t respectively. To further investigate this 
question and avoid potential systematic errors from the spectral 
method determination, the hardness ratios were derived (defined 
as the ratio of the excesses in different energy bands), using for 
this the energy (TeV) bands [0.2-0.35], [0.35-0.6] and [0.6-5.0]. 
For any combination, no differences were found beyond the 1 <x 
level between the increasing and decreasing parts. A similar ap- 
proach has been applied — when possible — for the rest of the 
flaring period. No clear dependence has been found within the 
statistical error limit of the determined indices, which is dis- 
tributed between 0.09 and 0.20. 

Finally, the persistence of the energy cut-off in the differ- 
ential energy spectrum along the flaring period has been exam- 
ined. For this purpose, runs were sorted again by increasing flux 
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Fig. 8. Integrated flux versus time for PKS 2155-304 on MJD 53944-53947 for four energy bands and with a 4-minute binning. 
From top to bottom: > 0.2 TeV, 0.2-0.35 TeV, 0.35-0.6 TeV and 0.6-5 TeV. These light curves are obtained using a power law 
spectral shape with an index of T = 3.37, also used to derive the flux extrapolation down to 0.2 TeV when the threshold is above that 
energy in the top panel (grey points). Because of the high dispersion of the energy threshold of the instrument (see Section|2] Fig.QTJi, 
and following the prescription described in l3.ll the integral flux has been determined for a time bin only if the corresponding energy 
threshold is lower than E m m. The fractional r.m.s. for the light curves are respectively, 0.86 ± 0.01 stat , 0.79 + 0.01 stat , 0.89 + 0.01 sta t 
and 1.01 ± 0.02 stat . The last plot shows the variation of the photon index determined in the 0.2-1 TeV range. See Section 14721 and 
Appendix lB.4l for details. 



and grouped into subsets containing at least an excess of 3000 y 
above 200 GeVQ. For the seven subsets found, the energy spec- 
trum has been determined in the 0.2-10 TeV energy range both 
for a simple power law and a power law with an exponential 
cut-off. This last hypothesis was found to be favoured systemati- 
cally at a level varying from 1.8 to 4.6 <x compared to the simple 
power-law and is always compatible with a cut-off in the 1-2 
TeV range. 

5. Light curve variability and correlation studies 

This section is devoted to the characterization of the temporal 
variability of PKS 2155-304, focusing on the flaring period ob- 
servations. The high number of y-rays available not only al- 
lowed minute-level time scale studies, such as those presented 
for MJD 53944 in lAharonian et al.l (2007a). but also to derive de- 



tailed light curves for three energy bands (Fig. [8]): 0.2-0.35 TeV, 
0.35-0.6 TeV and 0.6-5 TeV. 

The variability of the energy-dependent light curves of 
PKS 2155-304 is in the following quantified through their 
fractional r.m.s. F var defined in Eq. |4] (N andra et al. 19971 
lEdelson et al. 200% . In addition, possible time lags between 
light curves in two energy bands are investigated. 

5.1. Fractional r.m.s. F var 

All fluxes in the energy bands of Fig. [8] show a strong vari- 
ability which is quantified through their fractional r.m.s. F vm - 
(which depends on observation durations and their sampling). 
Measurement errors cr, elT on each of the N fluxes 0, of the light 
curve are taken into account in the definition of F var : 



5 To be significant, the determination of an energy cut-off needs a 
greater number of y than for a power law fit. 



(4) 
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where 5 2 is the variance 



1 N 

* 2 = ^2>-^ 



(5) 



;'=1 



and where <x 2 n . is the mean square error and (p is the mean flux. 

The energy-dependent variability F vm -{E) has been calcu- 
lated for the flaring period according to Eq.|4]in all three energy 
bands. The uncertainties on F vsir (E) have been estimated accord- 



ing to the parametrization derived by Vaughan et al. (2003), us- 
ing a Monte Carlo approach which accounts for the measure- 
ment errors on the simulated light curves. 

Fig. [9] shows the energy dependence of F vm - over the four 
nights for a sampling of 4 minutes where only fluxes with a sig- 
nificance of at least 2 standard deviations were considered. There 
is a clear energy-dependence of the variability (a null probability 
of ~ 10~ 16 ). The points in Fig.0 are fitted according to a power 
law showing that the variability follows F v . dr (E) oc £;0.i9±o.oi 
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Fig. 9. Fractional r.m.s. F var versus energy for the observation 
period MJD 53944-53947. The points are the mean value of the 
energy in the range represented by the horizontal bars. The line 
is the result of a power law fit where the errors on F wal and 
on the mean energy are taken into account, yielding F vm -(E) oc 

£0.19±0.01 



This energy dependence of F var is also perceptible within 
each individual night. In Table [6] the values of F var , the rela- 
tive mean flux, and the observation duration, are reported night 
by night for the flaring period. Because of the steeply falling 
spectra, the low energy events dominate in the light curves. This 
lack of statistics for high energy prevents to have a high fraction 
of points with a significance more than 2 standard deviation in 
light curves night by night for the three energy bands previously 
considered. On the other hand, the error contribution dominates 
and does not allow to estimate the F vm - in all these three energy 
bands. Therefore, only two energy bands were considered: low 
(0.2-0.5 TeV) and high (0.5-5 TeV). As can be seen in Table [6] 
also night by night the high-energy fluxes seem to be more vari- 
able than those at lower energies. 



Table 6. Mean Flux and the fractional r.m.s. F val night by night 
for MJD 53944-53947. The values refers to light curves with 4 
minute bins and respectively in three energy bands: >0.2 TeV, 
0.2-0.5 TeV, 0.5-5.0 TeV. Since a significant fraction (« 40%) 
of the points in the light curve of MJD 53945 in the energy band 
0.5-5.0TeV have a significance of less than 2 standard devia- 
tions, its F ml is not estimated. 



MJD 


Duration 


Energy 




p 




(min) 


(TeV1 


(10 _10 cnrV) 




53944 


88 












all 


15.44±0.87 


0.56+0.01 






0.2 - 0.5 


13.28±0.85 


0.55+0.01 






0.5-5.0 


1.94+0.24 


0.61+0.03 


53945 


244 












all 


2.40±0.41 


0.67+0.03 






0.2 - 0.5 


2.35±0.42 


0.64+0.03 






0.5-5.0 


0.34±0.12 




53946 


252 












all 


11.39+0.80 


0.35+0.01 






0.2 - 0.5 


10.02±0.79 


0.33+0.01 






0.5-5.0 


1.39±0.20 


0.43+0.02 


53947 


252 












all 


4.26+0.52 


0.22+0.02 






0.2 - 0.5 


4.02+0.52 


0.22+0.02 






0.5-5.0 


0.37+0.11 


0.13+0.09 



5.2. Doubling/halving timescale 

While F v - dr characterizes the mean variability of a source, the 



shortest doubling/halving time (Zhang et al. 1999) is an impor 



tant parameter in view of finding an upper limit on a possible 
physical shortest time scale of the blazar. 

If <D, represents the light curve flux at a time T,, for each 
pair of <D, one may calculate TV = \<S>AT/A<&\, where AT-Tj- 
T h A<t>=0 ; -cD, and O = (®y + 0,)/2. Two possible definitions 
of the doubling/halving are proposed by |Zhang et al.| (1999): the 
smallest doubling time of all data pairs in a light curve (T2), or 
the mean of the 5 smallest T'^ J (in the following indicated as 7^). 
One should keep in mind that, according to |Zhang et al. (1999), 
these quantities are ill defined and strongly depend on the length 
of the sampling intervals and on the signal-to-noise ratio in the 
observation. 

This quantity was calculated for the two nights with the high- 
est fluxes, MJD 53944 and MJD 53946, considering light curves 
with two different binnings (1 and 2 minutes). Bins with flux 
significances more than 2 <x and flux ratios with an uncertainty 
smaller than 30% were required to estimate the doubling time 
scale. The uncertainty on T2 was estimated by propagating the 
errors on the <£,, and a dispersion of the 5 smallest values was 
included in the error for T2. 

In Table UJ the values of T2 and T2 for the two nights are 
shown. The dependence with respect to the binning is clearly vis- 
ible for both observables. In this table, the last column shows that 
the fraction of pairs in the light curves which are kept in order to 
estimate the doubling times is on average ~45%. Moreover, dou- 
bling times T2 and T2 have been estimated for two sets of pairs 
in the light curves where A<t>=<l>j-<t>, is increasing or decreasing 
respectively. The values of the doubling time for the two cases 
are compatible within \,cr, therefore no significant asymmetry 
has been found. 

It should be noted that these values are strongly depen- 
dent on the time binning and on the experiment's sensitivity, so 
that the typical fastest doubling timescale should be conserva- 
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Table 7. Doubling/Halving times for the high intensity nights 
MJD 53944 and MJD 53946 estimated with two different sam- 
plings, using the two definitions explained in the text. The final 
column corresponds to the fraction of flux pairs kept to estimate 
the doubling times. 



MJD 


Bin size 


r 2 [min] 


T 2 [min] 


Fraction of pairs 


53944 


1 min 


1.65±0.38 


2.27±0.77 


0.53 


53944 


2 min 


2.20±0.60 


4.45±1.64 


0.62 


53946 


1 min 


1.61+0.45 


5.72±3.83 


0.25 


53946 


2 min 


4.55+1.19 


9.15±4.05 


0.38 



tively estimated as being less than ~ 2 min, which is compati- 
ble with the values reported in Ahar onian et al.l (2007a) and in 
lAlbert et aLl (2007). the latter concerning the blazar Mrk 501. 

5.3. Cross-correlation analysis as a function of energy 

Time lags between light curves at different energies can provide 
insights into acceleration, cooling and propagation effects of the 
radiative particles. 

The Discrete Correlation Function (DCF) as a function of the 
delay (IWhite & Peterson 19841 TEdelson & Krolik 1988b is used 
here to search for possible time lags between the energy-resolved 
light curves. The uncertainty on the DCF has been estimated us- 
ing simulations. For each delay, 10 light curves (in both en- 
ergy bands) have been generated within their errors, assuming 
a Gaussian probability distribution. A probability distribution 
function (PDF) of the correlation coefficients between the two 
energy bands has been estimated for each set of simulated light 
curves. The r.m.s. of these PDF are the errors related to the DCF 
at each delay. Fig.fTOlshows the DCF between the high and low 
energy bands for the four-night flaring period (with 4 minute 
bins) and for the second flaring night (with 2 minute bins). The 
gaps between each 28 minute run have been taken into account 
in the DCF estimation. 

The position of the maximum of the DCF has been estimated 
by a Gaussian fit, which shows no time lag between low and high 
energies for either the 4 or 2 minute binned light curves. This 
sets a limit of 14 ± 41 s from the observation of MJD 53946. A 
detailed study on the limit on the energy scale on which quantum 
gravity effects could become important, using the same data set, 
are reported in I Aharonian et al.l (2008 a). 

5.4. Excess r.m.s.-flux correlation 

Having defined the shortest variability time scales, the nature of 
the process which generates the fluctuations is investigated, us- 
ing another estimator: the excess r.m.s.. It is defined as the vari- 
ance of a light curve (Eq. [5) after subtracting the measurement 
error: 



0~x 



0"err)- 



(6) 



Fig. fTTI shows the correlation between the excess r.m.s. of 
the light curve and the flux, where the flux here considered are 
selected with an energy threshold of 200 Ge V. The excess vari- 
ance is estimated for 1- and 4-minute binned light curves, using 
20 consecutive flux points <t>, which are at least at the 2 cr signif- 
icance level (81% of the 1 minute binned sample). The correla- 
tion factors are r x = 0.60!° f 5 and r 4 = 0.87+° ->° for the 1 and 
4 minute binning, excluding an absence of correlation at the 2 cr 
and 4 cr levels respectively, implying that fluctuations in the flux 
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Fig. 10. Discrete Correlation Function (DCF) between the light 
curves in the energy ranges 0.2-0.5 TeV and 0.5-5 TeV and 
Gaussian fits around the peak. Full circles represent the DCF 
for MJD 53944-53947 4-minute light curve and the solid line is 
the Gaussian fit around the peak with mean value of 43 + 5 1 s. 
Crosses represent the DCF for MJD 53946 with a 2-minute light 
curve binning, and the dashed line in the Gaussian fit with a peak 
centred at 14 + 41 s. 



are probably proportional to the flux itself which is a characteris- 
tic of lognormal distributions (Aitchin son & Brown 19 63). This 
correlation has also been investigated extending the analysis to 
a statistically more significant data set including observations 
with a higher energy threshold in which the determination of the 
flux above 200 GeV requires an extrapolation (grey points in the 
top panel in Fig. [8]). In this case the correlations found are com- 
patible {r x = 0.78!°' \\ and r 4 = 0.93!° °^ for the 1 and 4 minute 
binning, respectively) and also exclude an absence of correlation 
with a higher significance (4 cr and 7 cr, respectively). 

Such a correlation has already been observed for 
X-rays in the Seyfert class A GN JEdelson et al. 20021 
|Vaughan, Fabian & Nandra 2003| |Vaughan et al 
|M c Hardy et al. 2004|> and - in 



2003, 
binaries 



and in X-ray 
( |Uttley & M c Hardy 200 1| |Uttley 2004| IGleissner et al. 2004K 
where it is considered as evidence for an underlying stochastic 
multiplicative process (Uttley, M c Hardy & Vaughan 2005 1, 
as opposed to an additive process. In additive processes, light 
curves are considered as the sum of individual flares "shots" 
contributing from several zones (multi-zone models) and the 
relevant variable which has a Gaussian distribution (namely 
Gaussian variable) is the flux. For multiplicative (or cascade) 
models the Gaussian variable is the logarithm of the flux. 
Hence, this first observation of a strong r.m.s.-flux correlation 
in the VHE domain fully confirms the log-normality of the flux 
distribution presented in Section l3~2l 



6. Characterization of the lognormal process during 
the flaring period 

This section investigates whether the variability of 
PKS 2155-304 in the flaring period can be described by a 
random stationary process, where, as shown in Section 15.41 
the Gaussian variable is the logarithm of the flux. In this 
case the variability can be characterized through its Power 
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will be further determined on the basis of parametric estimation 
and simulations. 

In the hypothesis where the process is stationary, i.e., the 
PSD is time-independent, a power law shape of the PSD was 
assumed, as for X-ray emitting blazars. The PSD was de- 
fined as depending on two paramenters and as follows: P(v) = 
K(v re f/v) a , where a is the variability spectral index and K de- 
notes the "power" (i.e., the variance density) at a reference fre- 
quency v le f . This latter was conventionally chosen to be 1CT 4 Hz, 
where the two parameters a and K are found to be decorrelated. 
Since a lognormal process is considered, P(v) is the density of 
variance of the Gaussian variable In <t>. The natural logarithm 
of the flux is conveniently used here, since its variance over a 
given frequency intervafl is close to the corresponding value of 
F^ a , at least for small fluctuations. For a given set of a and K, 
it is possible to simulate many long time series, and to modify 
them according to experimental effects, namely those of back- 
ground events and of flux measurement errors. Light curve seg- 
ments are further extracted from this simulation, with exactly the 
same time structure (observation and non-observation intervals) 
and the same sampling rates as those of real data. The distribu- 
tions of several observables obtained from simulations for dif- 
ferent values of a and K will be compared to those determined 
from observations, thus allowing these parameters to be deter- 
mined from a maximum-likelihood fit. 

The simulation characteristics will be described in Section [67TI 
Sections 16.21 16.31 and 16.41 will be devoted to the determi- 
nation of a and K by three methods, each of them based 
on an experimental result: the excess r.m.s.-flux correlation, 
the Kolmogorov first-order structure function ( IRutman 19781 
ISimonetti et al. 1985] ) and doubling-time measurements. 



Fig. 11. The excess r.m.s. cr xs vs mean flux <D for the observa- 
tion in MJD 53944-53947 (Full circles). The open circles are 
the additional points obtained when also including the extrapo- 
lated flux points - see text). Top: cr xs estimated with 20 minute 
time intervals and a 1 minute binned light curve. Bottom: cr xs 
estimated with 80 minute time intervals and a 4 minute binned 
light curve. The dotted lines are a linear fit to the points, where 
<x „ oc 0. 1 3 x <I> for the 1 minute binning and cr xs oc 0.3 x <t> for the 
4 minute binning. Fits to the open circles yield similar results. 



6.1. Simulation of realistic time-series 

For practical reasons, simulated values of In <£ were calculated 
from Fourier series, thus with a discrete PSD. The fundamental 
frequency vo = 1 /Tq, which corresponds to an elementary bin 
5v = vo in frequency, must be much lower than l/T if T is the 
duration of the observation. The ratio Tq/T was chosen to be 
of the order of 100, in such a way that the influence of a finite 
value of Tq on the average variance of a light curve of duration 
T would be less than about 2%. Taking Tq = 9 X 10 5 s, this 
condition is fulfilled for the following studies. With this approx- 
imation, the simulated flux logarithms are given by: 



Spectral Density (PSD) ( I van der Kris 19971) which indicates 
the density of variance as function of the frequency v. The 
PSD is an intrinsic indicator of the variability, usually rep- 
resented in large frequency intervals by power laws (oc v~ a ) 
and is often used to define different "states" of variable ob- 
jects (see e.g., IPaltani et al. 19971 and Zhang etal. 1999 for 
the PSD of PKS 2155-304 in the optical and X-rays). The 
PSD of the light curve of one single night (MJD 53944) was 
given in Aharonian et al. (2007a) between 10~ 4 and 10^ Hz, 
and was found to be compatible with a red noise process 
(a >2) with ~ 10 times more power as in archival X-ray 



data (Zhang etal. 1999), but with a similar index. This study 
implicitely assumed the y-ray flux to be the Gaussian variable. 
In the present paper, the PSD is determined using data from 
4 consecutive nights (MJD 53944-53947) and assuming a 
lognormal process. Since direct Fourier analysis is not well 
adapted to light curves extending over multiple days and 
affected by uneven sampling and uneven flux errors, the PSD 



In <D(f) = flo + ^ a„ cos(2njrvot + ip n ) 



(7) 



where N m - dx is chosen in such a way that 7o/W max is less than the 
time interval between consecutive measurements (i.e., the sam- 
pling interval). According to the definition of a Gaussian random 
process, the phases ip„ are uniformly distributed between and 
In and the Fourier coefficients a„ are normally distributed with 
mean and variances given by P(y) Sv/2 with v = ndv — n vq. 
From the long simulated time-series, light curve segments were 
extracted with the same durations as the periods of continuous 
data taking and with the same gaps between them. The simu- 
lated fluxes were further smeared according to measurement er- 
rors, according to the observing conditions (essentially zenith 
angle and background rate effects) in the corresponding data set. 



If a 2 is the variance of In O, F var = -^exp(cr 2 ) - 1. 
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6.2. Characterization of the lognormal process by the excess 
r.m.s.-flux relation 

For a fixed PSD, characterized by a set of parameters {a, K), 500 
light curves were simulated, reproducing the observing condi- 
tions of the flaring period (MJD 53944-53947), according to the 
procedure explained in Section l6"Tl 

For each set of simulated light curves, segments of 20 min- 
utes duration sampled every minute (and alternatively segments 
of 80 minutes duration sampled every 4 minutes) were extracted 
and, for each of them, the excess r.m.s. cr xs and the mean flux 
$ were calculated as explained in Section IB~4l For a large range 
of values of a and K, simulated light curves reproduce well the 
high level of correlation found in the measured light curves. On 
the other hand, the fractional variability F var and <I> are essen- 
tially uncorrected and will be used in the following. A likeli- 
hood function of a and K was obtained by comparing the simu- 
lated distributions of F val and O to the experimental ones. An ad- 
ditional observable which is sensitive to a and K is the fraction 
of those light curve segments for which F var cannot be calcu- 
lated because large measurement errors lead to a negative value 
for the excess variance. The comparison between the measured 
value of this fraction and those obtained from simulations is also 
taken into account in the likelihood function. The 95% confi- 
dence contours for the two parameters a and K obtained from 
the maximum likelihood method are shown in Fig. Q~2] for both 
kinds of light curve segments. The two selected domains in the 
{a, K) plane have a large overlap which restricts the values of a 
to the interval (1.9, 2.4). 



6.3. Characterization of the lognormal process by the 
structure function analysis 

Another method for determining a and K is based on 
Kolmogorov structure functions (SF). For a signal X(t), 
measured at N pairs of times separated by a delay r, 
[ti, tj + t} (i = 1, ...,N), the first-order structure function is de- 
fined as dSimonetti et al. 1985l >: 

1 N 

SF(T)^-Yj [X(ti) ~ X(ti + T)]2 (8) 

1=1 

1 N 

= ^2 [ln(1)(f '' )_ln0(r '' + r)]2 

In the present analysis, X{t) represents the variable whose PSD 
is being estimated, namely InO. The structure function is a 
powerful tool for studying aperiodic signals ( IRutman 19781 
ISimonetti et al. 1985b , such as the luminosity of blazars at var- 
ious wavelengths. Compared to the direct Fourier analysis, the 
SF has the advantage of being less affected by "windowing ef- 
fects" caused by large gaps between short observation periods in 
VHE observations. The first-order structure function is adapted 
to those PSDs whose variability spectral index is less than 3 
IRutman 19781 which is the case here, according to the results 
of the preceding section. 

Fig. [T3l shows the first-order SF estimated for the flaring pe- 
riod (circles) for r < 60 hr. At fixed r, the distribution of SF(r) 
expected for a given set of parameters {a, K) is obtained from 
500 simulated light curves. As an example, taking a — 2 and 
log 10 (/:/Hz _1 ) = 2.8, values of SF(t) are found to lie at 68% 
confidence level within the shaded region in Fig. [13] 



N 

^3.2 



g> 3.1 
O 



2.7 
2.6 













\ / 




' .-••■/:, 


\ / 










/ ■■ ; 














/ 




\ 






\ 











2.8 

a 



Fig. 12. 95% confidence domains for a and K at v re f = 10 _4 Hz 
obtained by a maximum-likelihood method based on the cr xs - 
flux correlation from 500 simulated light curves. The dashed 
contour refers to light curve segments of 20 minutes duration, 
sampled every minute. The solid contour refers to light curve 
segments of 80 minutes duration, sampled every 4 minutes. The 
dotted contour refers to the method based on the structure func- 
tion, as explained in Section l6~3l 
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Fig. 13. First order structure function SF for the observations 
carried out in the period MJD 53944-53947 (circles). The 
shaded area corresponds to the 68% confidence limits obtained 
from simulations for the lognormal stationary process charac- 
terized by a — 2 and log 10 (/T/Hz _1 ) = 2.8. The superimposed 
horizontal band represents the allowed range for the asymptotic 
value of the SF as obtained in Section|7] 

In the case of a power law PSD with index a, the SF aver- 
aged over an ensemble of light curves is expected show a vari- 
ation as t ([Kataoka et al. 2001 ). However, this property does 
not take into account the effect of measurement errors, nor of 
the limited sensitivity of Cherenkov telescopes at lower fluxes. 
For the present study, it was preferable to use the distributions 
of SF(t) obtained from realistic simulations including all exper- 
imental effects. Using such distributions expected for a given set 
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Table 8. Confidence interval at 68% c.l. for Tn and predicted 
by simulations for a— 2 and log 10 (.ST/ (Hz -1 )) = 2.8 for the two 
high-intensity nights MJD 53944 and MJD 53946, with two dif- 
ferent sampling intervals (1 and 2 minutes). 



MJD 


Bin size 


72 [min] 


Ti [min] 


53944 


1 min 


0.93-1.85 


1.60-2.60 


53944 


2 min 


3.01-4.28 


4.52-6.40 


53946 


1 min 


1.8-2.3 


1.96-2.41 


53946 


2 min 


5.3-9.1 


6.6-12.1 



of parameters {a, K], a likelihood function can be obtained from 
the experimental SF and further maximized with respect to these 
two parameters. Furthermore, the likelihood fit was restricted to 
values of t lower than 10 4 s, for which the expected fluctuations 
are not too large and are well-controlled. The 95% confidence 
region in the {a, K) plane thus obtained is indicated by the dot- 
ted line in Fig. [12] It is in very good agreement with those based 
on the excess r.m.s.-flux correlation and give the best values for 
a and K: 



a = 2.06 + 0.21 and log^/Hz -1 ) = 2.82 + 0.08 



(9) 



The variability index a at VHE energies is found to be re- 
markably close to those measured in the X-ray domain on 
PKS 2155-304, Mrk421, and Mrk501 dKataoka et al. 20011 

6.4. Characterization of the lognormal process by doubling 
times 

Simulations were also used to investigate if the estimator T2 
can be used to constrain the values of a and K. However, for 
a less than 2, no significant constraints on those parameters 
are obtained from the values of T2- For higher values of a, 
doubling times only provide loose confidence intervals on K 
which are compatible with the values reported above. This can 
be seen in Table [8] showing the 68% confidence intervals pre- 
dicted for T2 and To for a lognormal process with a-2 and 
log 10 (/T/(Hz -1 )) = 2.8, as obtained from simulation. Therefore, 
the variability of PKS 2155-304 during the flaring period can be 
consistently described by the lognormal random process whose 
PSD is characterized by the parameters given by Eq. [9] 



7. Limits on characteristic time of PKS 21 55-304 



In Section I5.2I the shortest variability time scale of 
PKS 2155-304 using estimators like doubling times have 
been estimated. This corresponds to exploring the high 
frequency behaviour of the PSD. In this section the lower 
(< 10~ 4 Hz) frequency part of the PSD will be considered, 
aiming to set a limit on the timescale above which the PSD, 
characterized in Section|6] starts to steepen to a > 2. A break in 
the PSD is expected to avoid infrared divergences and the time 
at which this break occurs can be considered as a characteristic 
time, from which physical mechanisms occurring in AGN could 
be inferred. 

Clearly the description of the source variability during the 
flaring period by a stationary lognormal random process is in 
good agreement with the flux distributions shown in Fig. [3] 
Considering the second Gaussian fit in the right panel of Fig. [3] 
the excess variance in the flaring regime reported in Table [9] al- 
though affected by a large error, is an estimator of the intrinsic 



Table 9. Variability estimators (definitions in Section [BTTT i rela- 
tive to In <t> both for the "quiescent" and flaring regime, as de- 
fined in Section l3~2l Experimental variance (line 1), error contri- 
bution to the variance (line 2), and excess variance (line 3). The 
latter is directly comparable to F^ ar . 







"Quiescent" regime 


Flarin; 


I regime 


°"exp 


of In® 


0.304 ± 0.040 


1.78 


± 0.27 


of In® 


0.169 ±0.053 


0.022 


± 0.005 




of In® 


0.135 ±0.067 


1.758 


± 0.273 



variance of the stationary process. It has been demonstrated that 

2 cr^ s represents the asymptotic value of the first-order structure 
function for large values of the delay t dSimonetti et al. 19851 1. 
On the other hand, as already mentioned, a PSD proportional to 
v~ a with a ss 2 cannot be extrapolated to arbitrary low frequen- 
cies; equivalently, the average structure function cannot rise as 
t"~' for arbitrarily long times. Therefore, by setting a 95% con- 
fidence interval on log 10 SF asympt = log 10 (2 <x^ s ) of [0.38,0.66] 
from Table [9] it is possible to evaluate a confidence interval on 
a timescale above which the average value of the structure func- 
tion cannot be described by a power law. Taking account of the 
uncertainties on a and K given by Eq. [9] leads to the 95% con- 
fidence interval for this characteristic time r c h a i of the blazar in 
the flaring regime: 

3 hours < r c ha r < 20 hours 

This is compatible with the behaviour of the experimental struc- 
ture function at times t > 10 4 s (Fig. [T3J, although the large 
fluctuations expected in this region do not allow a more accurate 
estimation. In the X-ray domain, characteristic times of the order 
of one day or less have been found for several blazars includ- 
ing PKS 2155-304 (IKataokaet al. 20011 ). The results presented 
here suggest a strong similarity between the PSDs for X-rays 
and VHE y-rays during flaring periods. 



8. Discussion and conclusions 

This data set, which exhibits unique features and results, is 
the outcome of a long-term monitoring program and dedi- 
cated, dense, observations. One of the main results here is 
the evidence for a VHE y-ray quiescent-state emission, where 
the variations in the flux are found to have a lognormal dis- 
tribution. The existence of such a state was postulated by 
IStecker & Salamon 19961 in order to explain the extragalac- 
tic y-ray background at 0.03-100 G eV detected by EGRET 
(IFichtel 19 96 Sree kumar et al. 1998b as coming from quiescent- 
state unresolved blazars. Such a background has not yet been de- 
tected in the VHE range, as it is technically difficult with the at- 
mospheric Cherenkov technique to find an isotropic extragalac- 
tic emission and even more to distinguish it from the cosmic -ray 
electron flux (Aharonia n et al. 2008bl >. In addition, the EBL at- 
tenuation limits the distance from which ~TeV y-rays can prop- 
agate to ~ 1 Gpc ([Aharonia n et al. 200 7b). As pointed out by 



Cheng, Zhang and Zhang (2000) emission mechanisms might 



be simpler to understand during quiescent states in blazars, 
and they are also the most likely state to be found observa- 
tionally. In the X-ray band, the existence of a steady underly- 
ing emission has also been invoked for two other VHE emit- 
ting blazars (Mrk 4 21, [Fo ssati et al. 200"0l and 1ES 1959+650, 
IGiebels et al. 20 02 ). Being able to separate, and detect, flaring 
and nonflaring states in VHE y-rays is hence important for such 
studies. 
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The observation of the spectacular outbursts of 
PKS 2155-304 in July 2006 represents one the most ex- 
treme examples of AGN variability in the TeV domain, and 
allows spectral and timing properties to be probed over two 
orders of magnitude in flux. 

Whereas for the flaring states with fluxes above a few 
10~ 10 cm' 2 s _1 a clear hardening of the spectrum with increasing 
flux is observed, familiar also from the blazars Mrk421 and Mrk 
501, for the quiescent state in contrast an indication of a soften- 
ing is noted. If confirmed, this is a new and intruiging observa- 
tion in the VHE regime of blazars. The blazar PKS 0208-512 
(of the FSRQ class) also shows such initial softening and sub- 
sequent hardening with flux in the MeV range, but no general 
trend could be found for y-ray blazars (Nandikot kur et al. 20071 ). 
In the framework of synchrotron self-Compton scenarios, VHE 
spectral softening with increasing flux can be associated with, 
for example, an increase in magnetic field intensity, emission 
region size, or the power law index of the underlying electron 
distribution, keeping all other parameters constant. A spectral 
hardening can equally be obtained by increasing the maximal 
Lorentz factor of the electron distribution or the Doppler factor 
(see e.g. Fig. 11.7 in lKataoka 19991 1. A better understanding of 
the mechanisms at play would require multi-wavelength obser- 
vations of similar time span and sampling density as the data set 
presented here. 

It is shown that the variability time scale f var of a few 
minutes are only upper limits for the intrinsic lowest char- 
acteristic time scale. Doppler factors of 6 > 100 of the 
emission region are derived by Ah aronian et al.l (2007a) using 
the ~ 10 9 M Q black hole (BH) Schwarzschild radius light 
crossing time as a limit, while Begelman et al. (2008) argue 



PKS 2155-304 the evidence of this log-normality has been 
found very recently in X-rays (Giebels & Degrange 2009 1 and as 
previously mentioned, X-ray binaries and Seyfert galaxies also 
show lognormal variability, which is thought to originate from 
the accretion disk (|M c Hardy et al. 2004| [Lyubarskii et al. 1997 



that such fast time scales cannot be linked to the size of 
the BH and must occur in regions of smaller scales, such 
as "needles" of matter moving faster than average within a 
larger jet (Ghisellini & Tavecchio 2008|, small components in 



the jet dominating at TeV energies ( Katarzyhski et al. 2008 
or jet "stratification" (Boutelie r, Henri & Petrucci 2008 



|Levinson (2007)| attributes the variability to dissipation in the jet 
coming from radiative deceleration of shells with high Lorentz 
factors. 

The flaring period allowed the study of light curves in sep- 
arated energy bands and the derivation of a power law depen- 
dence of F val with the energy (,F V ar K E~ 02 ). This depen- 
dence is comparable to that reported in IGiebels et al.l (2007), 
ILichti et al.l (2008), IMaraschi et all (2002), where F ml (E) oc 
E~ 02 between the optical and X-ray energy bands was found 
for Mrk 421 and PKS 2155-304, respectively. An increase with 
the energy of the flux variability has been found for Mrk 501 
( I Albert et al.l 2007) in VHE y-rays on timescales comparable to 
those observed here. 

The flaring period showed for the first time that the intrin- 
sic variability of PKS 2155-304 increases with the flux, which 
can itself be described by a lognormal process, indicating that 
the aperiodic variability of PKS 2155-304 could be produced 
by a multiplicative process. The flux in the "quiescent regime", 
which is on average 50 times lower than in the flaring period and 
has a 3 times lower F var , also follows a lognormal distribution, 
suggesting similarities between these two regimes. 

It has been possible to characterize a power spectral den- 
sity of the flaring period in the frequency range 10~ 4 -10 -2 Hz, 
resulting in a power law of index a = 2.06 + 0.21 valid for fre- 
quencies down to ~ 1/day. The description of the rapid vari- 
ability of a TeV blazar as a random stationary process must 
be taken into account by time-dependent blazar models. For 



|Arevalo and Uttley 2006 ), suggesting a connection between the 
disk and the jet. This variability behaviour should therefore be 
searched for in existing blazar light curves, independently of the 
observed wavelength. 
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Appendix A: Observations summary 

The journal of observations for the 2005-2007 is presented in 
TablerATTI 
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Table A.l. Summary of the 2005 to 2007 observations, where for each night MJD is the Modified Julian Date, 6 Z the mean 
observation zenith angle (degrees), T the total observation live-time (hours), « on the number of on-source events, n g the number 
of off-source events (from a region five times larger than for the on-source events). The final three columns are the corresponding 
excess, significance (given in units of standard deviations), and the significance per square root of the live-time. 



MJD 


0t 


T 


n on 


"off 


Excess 




o-/Vr 


53618 


16.9 


0.87 


860 


3,202 


219.6 


7.5 


8.0 


53637 


16.6 


1.09 


788 


2,673 


253.4 


9.2 


8.8 


53638 


14.7 


2.18 


1,694 


6,029 


488.2 


12.0 


8.1 


53665 


22.6 


0.87 


618 


2,487 


120.6 


4.7 


5.1 


53666 


26.3 


1.31 


857 


3,406 


175.8 


5.9 


5.1 


53668 


23.3 


1.30 


926 


3,793 


167.4 


5.3 


4.7 


53669 


19.6 


0.86 


1,027 


2,939 


439.2 


14.7 


15.8 


53705 


55.6 


0.88 


512 


2,542 


3.6 


0.1 


0.2 


53916 


13.7 


0.88 


993 


3,317 


329.6 


10.7 


11.5 


53917 


11.8 


0.88 


933 


3,163 


300.4 


10.1 


10.7 


53918 


10.2 


1.32 


1,491 


4,596 


571.8 


15.5 


13.5 


53919 


10.9 


1.32 


1,477 


4,638 


549.4 


14.9 


13.0 


53941 


14.4 


1.31 


2,445 


4,844 


1,476.2 


35.0 


30.6 


53942 


13.7 


1.76 


2,453 


5,766 


1,299.8 


29.5 


22.3 


53943 


9.8 


1.33 


1,142 


3,627 


416.6 


12.8 


11.1 


53944 


13.2 


1.33 


12,762 


3,563 


12,049.4 


172.9 


149.7 


53945 


23.9 


5.23 


8,037 


16,352 


4766.6 


62.0 


27.1 


53946 


27.7 


6.61 


35,874 


19,881 


31,897.8 


251.3 


97.7 


53947 


25.1 


5.89 


17,158 


17,006 


13,756.8 


142.6 


58.8 


53948 


27.7 


2.75 


5,366 


7,957 


3,774.6 


64.6 


38.9 


53950 


26.6 


3.51 


5,108 


11,955 


2,717.0 


42.8 


22.9 


53951 


28.4 


2.51 


3,275 


8,421 


1,590.8 


30.6 


19.3 


53952 


35.8 


1.76 


1,786 


5,395 


707.0 


17.7 


13.3 


53953 


44.1 


0.89 


670 


2,285 


213.0 


8.4 


8.9 


53962 


27.6 


0.89 


534 


2,088 


116.4 


4.9 


5.3 


53963 


19.4 


1.75 


1,613 


6,145 


384.0 


9.5 


7.1 


53964 


10.3 


1.49 


1,057 


4,146 


227.8 


6.9 


5.6 


53965 


15.7 


1.57 


1,584 


5,662 


451.6 


11.4 


9.1 


53966 


18.6 


0.88 


719 


2,844 


150.2 


5.5 


5.8 


53967 


24.3 


0.86 


481 


1,801 


120.8 


5.5 


5.9 


53968 


19.1 


0.86 


479 


1,974 


84.2 


3.7 


4.0 


53969 


21.2 


1.29 


1,738 


4,368 


864.4 


23.0 


20.2 


53970 


20.9 


0.88 


690 


2,759 


138.2 


5.1 


5.5 


53971 


19.8 


1.32 


1,449 


4,313 


586.4 


16.3 


14.2 


53972 


16.5 


1.30 


683 


2,499 


183.2 


7.0 


6.1 


53973 


14.5 


1.32 


1,157 


4,311 


294.8 


8.6 


7.5 


53974 


16.2 


1.76 


1,504 


5,925 


319.0 


8.1 


6.1 


53975 


15.8 


0.88 


804 


3,059 


192.2 


6.7 


7.2 


53976 


13.8 


0.89 


727 


2,544 


218.2 


8.2 


8.7 


53977 


12.0 


0.88 


832 


2,745 


283.0 


10.1 


10.8 


53978 


13.3 


0.89 


687 


2,317 


223.6 


8.7 


9.3 


53995 


21.4 


1.75 


1,712 


5,989 


514.2 


12.6 


9.5 


53996 


25.3 


1.33 


680 


2,834 


113.2 


4.2 


3.6 


53997 


26.3 


1.32 


1,113 


4,382 


236.6 


7.0 


6.1 


53998 


26.1 


1.32 


1,247 


4,082 


430.6 


12.6 


11.0 


53999 


21.4 


1.32 


1,107 


4,262 


254.6 


7.5 


6.6 


54264 


8.5 


0.13 


109 


449 


19.2 


1.8 


5.0 


54265 


8.5 


1.05 


920 


3,513 


217.4 


7.1 


6.9 


54266 


9.3 


1.39 


1,129 


4,553 


218.4 


6.3 


5.4 


54267 


9.2 


1.32 


1,095 


4,176 


259.8 


7.8 


6.7 


54268 


10.2 


0.32 


261 


1,040 


53.0 


3.2 


5.6 


54269 


9.1 


0.89 


908 


2,491 


409.8 


14.7 


15.6 


54270 


10.8 


0.44 


565 


1,266 


311.8 


14.9 


22.4 


54271 


7.6 


0.36 


308 


983 


111.4 


6.6 


11.0 


54294 


7.7 


0.44 


447 


1,337 


179.6 


9.0 


13.5 


54296 


7.2 


0.44 


350 


1,289 


92.2 


4.9 


7.4 


54297 


9.8 


0.44 


344 


1,265 


91.0 


4.9 


7.4 


54299 


7.6 


0.44 


347 


1,232 


100.6 


5.5 


8.2 


54300 


7.3 


0.44 


343 


1,245 


94.0 


5.1 


7.6 


54302 


7.7 


0.44 


358 


1,234 


111.2 


6.0 


9.0 


54304 


7.9 


0.88 


680 


2,765 


127.0 


4.7 


5.0 


54319 


8.9 


0.88 


692 


2,650 


162.0 


6.1 


6.5 


54320 


8.0 


0.89 


553 


2,199 


113.2 


4.7 


5.0 


54329 


11.7 


0.44 


297 


1,258 


45.4 


2.5 


3.8 


54332 


7.0 


0.16 


100 


391 


21.8 


2.1 


5.3 


54375 


9.4 


1.27 


811 


3,124 


186.2 


6.4 


5.7 


54376 


7.9 


0.68 


395 


1,605 


74.0 


3.6 


4.4 
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Appendix B: Spectral variability 

The numerical information associated with Fig. [7j is given in 
Tables lBTK left panel, grey points), lB~2l (left panel, black points) 
and lB.3l (right panel). In addition, numerical information associ- 
ated with Fig.[8]is given in Table IB .41 

Table B.l. Integral flux (KT 11 cirrV 1 ) in the 0.2-1 TeV energy 
range versus photon index corresponding to grey points in the 
left panel of Fig. [JJ Errors are statistical. See Section 14.11 for 
more details. 





Index T 


2.36 ±0.13 


3.345 ± 0.20 


3.92 ±0.17 


3.64 ±0.16 


5.33 ± 0.22 


3.46 ±0.13 


8.29 ± 0.30 


3.64 ±0.10 


13.82 ±0.82 


3.82 ±0.17 



Table B.2. Integral flux (10" n cirrV 1 ) in the 0.2-1 TeV energy 
range versus photon index corresponding to black points in the 
left panel of Fig. [JJ Errors are statistical. See Section 14.11 for 
more details. 



<D 




Index T 


8.4 ± 


0.3 


3.74 ± 


0.11 


16.9 ± 


0.5 


3.82 ± 


0.10 


24.5 ± 


0.7 


3.78 ± 


0.08 


37.4 ± 


1.0 


3.77 ± 


0.08 


39.7 ± 


1.1 


3.76 ± 


0.08 


46.4 ± 


1.1 


3.66 ± 


0.08 


53.5 ± 


1.3 


3.57 ± 


0.07 


78.6 ± 


1.9 


3.44 ± 


0.06 


91.8 ± 


1.9 


3.33 ± 


0.06 


101.6 d 


:2.8 


3.30 ± 


0.07 


111.7 d 


:3.0 


3.33 ± 


0.07 


154.1 d 


:3.5 


3.28 ± 


0.06 


173.1 d 


:3.8 


3.16 ± 


0.06 


198.5 d 


:3.8 


3.28 ± 


0.05 


210.9 ±3.9 


3.14 ± 


0.05 



Table B.4. MJD, integral flux (10 _n cirrV 1 ) in the 0.2- 
1 TeV energy range, and photon index corresponding to the en- 
tries of Fig [8] Only points associated with an energy threshold 
lower than 200 GeV are considered. Errors are statistical. See 
Section l4~2l for more details. 





MID 






<D 




Index T 


53944 


.02742 ± 





,00277 


188.6 ± 


30.6 


3.22 


±0.09 


53944 


.03298 ± 





,00277 


184.1 ± 


30.9 


3.28 


±0.09 


53944 


.03854 ± 





,00277 


191.7 ± 


32.6 


3.45 


±0.09 


53944 


.04409 ± 


0, 


,00277 


252.4 ± 


40.5 


3.19 


±0.09 


53944 


.04965 ± 





,00277 


237.5 ± 


33.8 


3.16 


±0.08 


53944 


,05520 ± 





,00277 


212.8 ± 


30.8 


3.04 


±0.08 


53944 


06076 ± 





,00277 


190.9 ± 


30.0 


3.09 


±0.09 


53944 


06909 ± 





,00555 


99.5 ± 


17.5 


3.18 


±0.10 


53944 


.98298 ± 





,05277 


9.2 ± 


3.1 


3.89 


±0.18 


53945 


.04965 ± 





,01388 


34.2 ± 


8.3 


3.83 


±0.13 


53945 


.07604 ± 





,01250 


42.6 ± 


10.2 


3.92 


±0.13 


53945 


.93020 ± 





,00277 


206.9 ± 


34.9 


3.2: 


t0.10 


53945 


.93715 ± 





00416 


190.9 ± 


28.6 


3.15 


±0.09 


53945 


.94409 ± 





,00277 


171.0 ± 


30.5 


3.23 


±0.10 


53945 


.94965 ± 





,00277 


161.2 ± 


28.2 


3.06 


±0.10 


53945 


.95659 ± 





00416 


173.9 ± 


28.9 


3.35 


±0.09 


53945 


.96354 ± 





,00277 


179.1 ± 


29.4 


3.11 


±0.10 


53945 


,97048 ± 





,00416 


127.8 ± 


22.6 


3.36 


±0.10 


53945 


,98159 ± 





00694 


91.9 ± 


16.5 


3.42 


±0.10 


53945 


,99687 ± 





,00833 


101.6 ± 


16.9 


3.12 


±0.10 


53946 


,00937 ± 





00416 


104.4 ± 


19.7 


3.28 


±0.11 


53946 


,01909 ± 





,00555 


92.2 ± 


17.1 


3.37 


±0.10 


53946 


,03020 ± 





,00555 


79.5 ± 


15.9 


3.42 


±0.11 


53946 


.03992 ± 





00416 


105.6 ± 


19.3 


3.26 


±0.10 


53946 


,04965 ± 





,00555 


114.7 ± 


20.9 


3.45 


±0.10 


53946 


,05937 ± 





00416 


110.2 ± 


19.6 


3.17 


±0.10 


53946 


06909 ± 





,00555 


97.2 ± 


19.9 


3.54 


±0.11 


53946 


,08298 ± 





,00833 


75.1 ± 


12.8 


3.24 


±0.10 


53946 


,09826 ± 





00694 


75.1 ± 


15.5 


3.55 


±0.12 


53946 


.93437 ± 





,00972 


58.1 ± 


13.5 


3.92 


±0.13 


53946 


.95381 ± 





,00972 


50.7 ± 


11.2 


3.6: 


t0.12 


53946 


.97604 ± 





,01250 


39.7 ± 


8.7 


3.62 


±0.12 


53947 


.00242 ± 





,01388 


34.0 ± 


7.5 


3.62 


±0.12 


53947 


.02742 ± 





01111 


40.4 ± 


8.8 


3.56 


±0.12 


53947 


,04965 ± 





01111 


43.5 ± 


10.3 


3.82 


±0.13 


53947 


.07742 ± 





01666 


45.8 ± 


10.1 


3.75 


±0.12 



Table B.3. Integral flux (10" n cirrV 1 ) in the 0.2-1 TeV energy 
range versus photon index corresponding to the right panel of 
Fig. [7] Errors are statistical. See Section l4~T1 for more details. 



O Index T 



4.2 ±0.1 


3 


.52 ± 


0.07 


17.3 ±0.5 


3 


.89 ± 


0.09 


38.6 ± 1.1 


3 


.80 ± 


0.08 


43.6 ± 1.1 


3 


.60 ± 


0.07 


51.5 ± 1.3 


3 


.53 ± 


0.07 


67.2 ± 1.9 


3 


,64 ± 


0.07 


86.1 ± 1.9 


3 


.38 ± 


0.06 


97.0 ± 1.9 


3 


,30 ± 


0.05 


111.7 ± 3.0 


3 


,33 ± 


0.07 


154.1 ±3.5 


3 


,28 ± 


0.06 


173.1 ±3.8 


3 


,16± 


0.06 


198.5 ±3.8 


3 


,28 ± 


0.05 


210.9 ±3.9 


3 


,14 ± 


0.04 



